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Results for the late-time regime of phase ordering in three dimensions are reported, based on 
numerical integration of the time-dependent Ginzburg-Landau equation with nonconserved order 
parameter at zero temperature. For very large systems (700 ) at late times, t > 150, the char- 
acteristic length grows as a power law, R(t) ~ t n , with the measured n in agreement with the 
theoretically expected result n = 1/2 to within statistical errors. In this time regime R(t) is found 
to be in excellent agreement with the analytical result of Ohta, Jasnow, and Kawasaki [Phys. Rev. 
Lett. 49, 1223 (1982)]. At early times, good agreement is found between the simulations and the 
linearized theory with corrections due to the lattice anisotropy. 
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I. INTRODUCTION 



Phase ordering and phase separation of materials, following a rapid change in an intensive variable from a region 
of the phase diagram where the system is uniform to one in which two or more phases coexist, are among the oldest 
and most common methods of materials processing. A typical example is the temperature quenching performed by 
blacksmiths since antiquity, in which hot metal is suddenly cooled by immersion in water. In fact, the metallurgical 
term "quenching" has become common in the literature on the dynamics of phase transformations. Modern examples 
of the use of phase ordering as a processing technique include precipitation strengthening in metals B and fabrication 

^ ■ of glasses @. 

As the domains of different phases evolve and grow after the quench, the dynamic scaling hypothesis states that 

rf\ • their behavior over a large range of length scales can be described in terms of a single, time-dependent characteristic 

(T) ' length, R(t). For many phase-ordering processes, this characteristic length behaves as a power law for asymptotically 

CN| . late times, 

t> 

O; R(t)~t n , (l) 



where the growth exponent n depends on the dynamic universality class MM- The simplest of these universality 
classes is comprised of systems with only local relaxational dynamics and a nonconserved scalar order parameter, 
known as "Model A" in the classification scheme of Hohenberg and Halperin |J . At late times the order parameter 
takes distinct values in the two phases, which without loss of generality can be taken as ±1. The two phases are 
separated by a sharp interface, where the order parameter is near zero, and the local interface velocity is proportional 
to the local mean curvature, H(r,t). In two dimensions this is simply the inverse of the local radius of curvature, 
i?(r, t), while in three dimensions it is the arithmetic mean, H (r, t) — [1/Ri(r, t) + l/R-zir, t)] /2, where R\ and i?2 are 
the two principal curvatures. The global characteristic length, R(t), can be identified as proportional to the inverse 
J> \ of the average of \H(r, t)\ over the whole interface. It thus obeys the asymptotic equation of motion, 



R(t) ~ 1/R(t) , (2) 

which yields n = 1/2, independent of the spatial dimension. This result was shown early on by Lifshitz fa], Chan 
||, and Allen and Cahn 0, and it is often referred to as Lifschitz-Allen-Cahn dynamics. Physical realizations of this 
universality class include phase ordering in anisotropic magnets B, alloys such as CU3AU M and FeaAl M, liquid 
crystals |1(|[JT| , and adsorbate systems |1J] . 

Since experimental complications due to other effects, such as strain fields and hydrodynamics, usually cannot 
be completely excluded, it is desirable to obtain numerical verification in a cleanly defined three-dimensional model 
system. Even in two dimensions direct numerical verification of the asymptotic t 1 ' 2 growth through a direct estimate 
of R(t) is uncommon; examples are Refs. |13-[l5|. A more common practice is to show consiste ncy of numerical results 



with the asymptotic growth, e.g. through Monte-Carlo Renormalization Group techniques |16|17|| or scaling of the 



correlation function |13] or structure factor |l9[ . However, to our knowledge only experimental verifications have so 
far been reported in three dimensions [g|,y_l| . Until now, numerical verification has been prevented by the very large 
systems and long simulation times needed to observe the asymptotic scaling over a sufficient time interval to provide 
accurate measurements of n. In this paper we present such unequivocal confirmation of t 1 ' 2 growth at late times in 
three-dimensional numerical simulations of very large systems. 

The remainder of this paper is organized as follows. In Sec. Owe introduce the numerical model and discuss the 
numerical method used to integrate its time evolution. In Sec. |EH| we discuss the time evolution at early times. In 
Sec. |V| we present the main numerical results of this paper: the growth of the characteristic length at late and 
intermediate times. A summary and conclusions are presented in Sec. M. 

II. MODEL AND NUMERICAL METHOD 

A generic model for the nonconserved dynamics of Model A is given by the time-dependent Ginzburg-Landau 
(TDGL) equation, 



dj>{v,t) _ 6FW(r,t)] 
dt Sip(r,t) 



CM), (3) 



where the functional derivative corresponds to the deterministic relaxation associated with the free-energy functional 
^[^(rjt)], and C( r >^) 1S a stochastic process that represents thermal fluctuations. For the local part of the free-energy 
functional we choose the Ginzburg-Landau- Wilson free energy |J] 



F[i>(r,t)]=[dr 



-^> 2 M) + ^ 4 (M) + ||v^M)| 2 



(4) 



which has minima at ip — ±1, corresponding to the two degenerate uniform phases. The problem can be cast into 
this dimensionless form without loss of generality papO] . The nonequilibrium process associated with a system that 
is quenched to a temperature far below the critical temperature is controlled by a zero-temperature fixed point [|l|. 
Thus, the stochastic part of Eq. (H) can be ignored, and the equation of motion becomes 

^A = (l + cV 2 )^r,t)-^(r 1 t). (5) 

The numerical integration of Eq. (jq) with c = 3/2 pl[ was performed using a finite-difference approach on cubic lattices 
with periodic boundary conditions and 100 3 , 300^500 3 , and 700 3 points. Results for each system size were averaged 
over 5 integrations from different initial conditions, except for the 300 3 lattice, for which results were averaged over 10 
runs. The initial condition consisted of a random value of i\) at each lattice point, with values chosen from a uniform 
distribution on [—tJiq^t/jq]. Unless otherwise noted, ipo=0.l for results presented here. From this initial configuration, 
the system was integrated using a first-order Euler scheme with Ai=0.01. For early times, t < 10, we also tried At 
reduced by a factor of 5, which changed R(t) by only about 2%. Approximate isotropy was ensured by using a 19-point 
discretization of the Laplacian, analogous to the 9-point discretization commonly used in two dimensions |23|23J . 

The computational resources required for this numerical integration are large. The storage required for one array 
of 700 3 lattice sites is more than 2.5 Gigabytes, and two such arrays are required by the integration algorithm. 
Integration of the 700 3 lattice over 1500 time units took 84 hours on 15 four-processor nodes (5040 CPU hours) on 
an IBM SP3 supercomputer. 

III. EARLY-TIME BEHAVIOR 

Immediately after the quench, the local order parameter ip (r, t) is randomly distributed with values centered around 
0. The initial response of the system is to form small regions in local equilibrium, dominated by values near +1 or —1. 
This process is essentially completed by t ss 10, as illustrated in Fig. [j], which shows the time evolution of y/(?p 2 (r, £)) 
for early times. For the initial condition used here (ip 2 (r, 0)) = ip 2 /3 immediately after the quench, and it approaches 
unity at late times as the regions in local equilibrium come to dominate the system. 

The initial relaxation away from the uncorrelated random state, towards a state dominated by regions in which the 
order parameter everywhere has the same sign, is well described by the linearized version of Eq. (pt) (corresponding 



to a Cahn-Hilliard equation for nonconserved order parameter |24|,|25|] ) . The Fourier representation of the solution of 
this linear dynamical equation is 

ip(k,t) =^(k,0)exp[(l-c£; 2 )£] , (6) 

where ip (k, t) is the spatial Fourier transform of the order-parameter field. The progress of the phase ordering at 
early times can be quantified by (ip 2 (r,t)}, where () represents averaging over space. By integrating "0(k, t)-0(— k, i) 
over the k-space region associated with the finite system, the spatial average can be evaluated as 

(V' 2 (r ) *)) = ^^-/dkexp(2t[l- C fe 2 (k)]) ! (7) 

where d is the spatial dimension, and the integration limits are [— ir, it] in each Cartesian coordinate. In deriving 
this result, we used the fact that the uncorrelated initial condition used here gives V»(k, 0)r/>(— k, 0) = (ip 2 (r,0)), 
independent of k. In Eq. (fil) the expression k 2 (k) takes into account the anisotropy of k 2 that results from the 
numerical implementation of the Laplacian. When no anisotropy exists, the integral in Eq. (f7|) can be evaluated 
analytically to yield 



<<A 2 (r,i)}HV> 2 (r,0))e 



erf(7rv/2ct) 



2%/2ttc£ 



(8) 



This result is shown as dashed curves in Fig. |lj. For the Laplacian used in the simulations presented here, the 
anisotropy in the squared magnitude of the wave vector is 



fc 2 (k)-4-l 



cos k x cos 2 — + cos k v cos 2 — + cos k z cos 2 — 
2 y 2 2 



(9) 



Using this expression for /c 2 (k), Eq. (R) was evaluated numerically using midpoint integration at 10 6 uniformly 
distributed points. The results are presented as the solid curves in Fig. [y, where good agreement is seen between the 
simulations and the linear theory for t < 5. The initial decrease is due to the diffusional decay of high-wavevector 
modes with fc 2 (k) > cr l . During this brief period one can define a microscopic diffusion length which increases 
with time as t 1 / 2 p3]. The subsequent rapid increase is caused by the exponential growth of the modes at smaller 
wavevectors. Similar time dependence for (ip 2 (r,t)} has been observed previously in two-dimensional simulations 
G3. As t/>(r, t) approaches unity, the neglected cubic term in the TDGL equation becomes important, and the linear 
approximation breaks down as the order parameter saturates to its degenerate equilibrium values inside the domains. 
To quantify the effects of the initial conditions on the early-time behavior, results for one simulation with ipo=l on 
a 300 3 system are also shown in Fig. [j]. Even for this large value of ipo, the linear theory works reasonably well at the 
earliest times. 

IV. LATE AND INTERMEDIATE TIMES 

The scaling ansatz associated with Eq. (|lj) is not valid until a clear separation has been achieved between large- 
scale fluctuations representing domains in which the order parameter takes values near its two degenerate equilibrium 
values, and microscopic fluctuations of the local order parameter about these values within the domains J9|,|o] p8| . 
While the early-time growth is completed by i«10, the separation of length scales is not complete until a significantly 
later time, as is now shown. 

For the large simulations presented here, it was necessary to use a computationally efficient estimate of the charac- 
teristic length R(t). This was done by identifying R(t) as proportional to the inverse of the interfacial area per unit 
volume [p9|,p0| . The interface area was measured by counting the number of the nearest-neighbor lattice-site pairs 
having values of the local order parameter with opposite signs. As a result, 

R=h/(3N)J2@[-^(r i )^(r j )]\ -1. (10) 

Here N is the number of lattice sites, the sum is over all nearest- neighbor pairs, and the Heaviside function Q[x] =0 
for x < and 1 otherwise. This definition of R{t) approximates the inverse r-derivative of the normalized two-point 



correlation function, C(r,t) = C(r/R(t)), in the small-r limit. Details on the derivation of Eq. ( JIOJ ) are given in the 
Appendix. Direct comparison with the much more computationally intensive C(r,t) for a 315 3 system confirms this 
equality to within 4% for t > 500. Corrections to R, related to unequal volumes of the degenerate phases p9| , p0[ , 
which become important at very late times, do not affect the results presented here and have been neglected. The 
measured values of R are shown in Fig. vs time on a log-log scale. The characteristic length clearly does not obey a 
single power law for all times. It is only at the latest times, t > 150, that the asymptotic power-law regime is reached. 
Least-squares fitting of a power law to the 700 3 data gives an exponent 0.511 ± 0.01 for 150 < t < 1500. This result 
is consistent with the expected value n = 1/2. 

A more sensitive test for true power-law behavior can be made by measuring the instantaneous, effective growth 
exponent as a function of time. Here this is accomplished by estimating the derivative d\nR/d\rit using 3-point 
central differencing. The results are presented in Fig. 0. During the early-time regime of near-exponential growth 
of ^ (tp 2 (r , t)) , the effective exponent for R(t) falls steeply from near unity at very early times to near 0.40 around 
t w 10. For t > 20 the effective exponent rises again for the systems larger than 100 3 . This steady increase of the 
effective exponent in the intermediate-time regime indicates that here, too, the dynamics are not properly described 
by a power law. The only system for which the exponent does become constant is the 700 3 lattice for t > 150. For 
these late times the mean exponent is 0.519 ± 0.01. The error on all estimates of the exponent reported here is the 
standard deviation of these data. Given the trends with system size observed in Figs. |2j and |3J, finite-size effects are 
most likely affecting the late-time behavior in the other system sizes considered here. 

A final test for n — 1/2 at late times is presented in Fig. [|, where R is plotted against y/t. The straight line is 
the least-squares fit of the 700 3 data for t > 150; it has a correlation coefficient of 0.999976. The effect of system 
size on the growth shows a clear trend, with progressively smaller systems departing from the common behavior at 
progressively earlier times. 

It is informative to compare the present results for R(t) to the theoretical prediction of Ohta, Jasnow, and Kawasaki 
(OJK) pll. Defined such that the slope of the normalized two-point correlation function C(r,t) is —1/R(i) in the 
small-r limit, the OJK result is 



RojK(t) = ^4c^-t. (11) 

In Fig. g, the OJK result appears as the dashed line, and the agreement with the simulation data at late times is 
excellent. The OJK theory is similarly successful for the two-dimensional analog of the model presented here, where 
both the agreement between Eqs. ( JIOJ ) and (|ll|) and the resulting scaled form of the structure factor (the Fourier 
transform of the two-point correlation function) are excellent |lij . 

An alternative method to obtain the specific interface area in this system, is to consider the quantity A(t) = 
1 — ^/(ifj 2 (r, £)), which corresponds to the interface area multiplied with an average interface thickness. Once the 
interface thickness has converged to a time- independent value, A~ l (t) oc R(t) [p3[ . For times later than approximately 
20, A^ 1 increases as a power law with t, as shown in Fig. pi Least-squares fits to the simulation data for 150 < t < 1500 
result in estimated exponents of 0.487±0.01, 0.480±0.01, and 0.512±0.01 for the 300 3 , 500 3 , and 700 3 , respectively. An 
average over the same interval of effective exponents for A _1 (£), obtained by 3-point differencing in a way analogous to 
those discussed above for R(t), give an estimate of 0.511 ±0.01 for the 700 3 system. These 700 3 results, in particular, 
are thus consistent with n=l/2. 

V. SUMMARY AND CONCLUSIONS 

From the numerical data obtained in this study, we find that the time evolution of the three-dimensional Model 
A system defined by Eq. (H) can be divided into four time regimes. The time regimes and the evolution behavior in 
each are summarized as follows. 

1. For very early times, before the order parameter has saturated to near its two degenerate equilibrium values 
inside distinct domains, the time evolution is well described by the linearized version of Eq. (||), as seen in Fig. [I]. 
During this stage, and until t sa 10, the effective exponent of R(t) falls steeply from near unity to near 0.45, as 
seen in Fig. |3[ 

2. In the intermediate-time regime, for 10 < t < 150, both R(t) and A~ l {t) as they are defined here become 
reasonable measures of a characteristic length, and least-squares fits to log-log plots, such as Fig. || and Fig. ||, 
in this time regime yield apparent exponent estimates near 0.45. However, inspection of Fig. shows that the 
effective exponent is not constant in this regime: for the systems larger than 100 3 it increases steadily back 



towards the vicinity of 0.5. Thus, the growth of the characteristic length in this intermediate-time regime is also 
clearly not well described by a simple power law. 

3. For late times, t > 150, and for the largest system studied, 700 3 , the effective exponent for R(t) [and also for 
A _1 (i)] levels off to fluctuate near 0.5. A least-squares fit to R(t) (Fig. ||) over a full decade, 150 < t < 1500, 
yields an estimate of 0.511 ± 0.01, while an average over the effective exponents (Fig. ||) in the same interval 
yields 0.519 ± 0.01. The corresponding estimates for A~ l (t) are 0.512 ± 0.01 and 0.511 ± 0.01, respectively. 
These estimates are all consistent with the theoretical expectation of ro=l/2. 

4. For some t > 1500 for the 700 3 system, and indeed for much smaller times for the smaller systems, finite-size 
effects made observation of t 1 ' 2 growth impossible. In this very-late-time regime, which is pushed out to later 
times for larger systems, R(t) becomes on the order of the system size, and the order parameter selects one or 
the other of its two degenerate values. In Fig. ||, and even more clearly in Fig. [|, it can be seen how R(t) for 
progressively smaller systems deviates from the asymptotic behavior at progressively earlier times. 

We emphasize that, although we see four different time regimes in the domain growth, it is only in the late-stage 
Lifshitz-Allen-Cahn regime (150 < t < 1500 for the 700 3 system) that true power-law growth is observed. While naive 
least-squares fits to Fig. |j indeed yield apparent exponents at early and intermediate times, similar to those recently 
published by Fialkowski et at Jfq], those regimes are not properly described by simple power laws. We also observe 
from our data that, in disagreement with what is claimed in Rcf. |15| , the t 1 ' 2 growth in general is observed at times 
before the final deviation of the average order parameter from zero. This final loss of symmetry is a finite-size effect, 
and it occurs only in the very-late-time regime. 

In conclusion, we have presented the first unequivocal results confirming t 1 ' 2 domain growth for integration of a 
three-dimensional numerical instance of Model A, describing phase ordering in a system with nonconserved order 
parameter. In this late-time regime of t 1 ' 2 growth we found excellent agreement between the observed characteristic 
length and the analytic result of Ohta, Jasnow, and Kawasaki [BlJ. In order to obtain these solid numerical results, 
very long simulations of very large systems were necessary. 
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APPENDIX 



In this Appendix we sketch the derivation of Eq. (nfl) for the characteristic length R(t) in terms of the total number 
of bonds between positive and negative nearest-neighbor ip(vi). 

In a (i-dimensional system with infinitely thin, randomly oriented interfaces, the inverse r-derivative of the normal- 
ized order-parameter correlation function C(r, t) is p9| , |30] | 

R(t) = (l-^) 2 )V/{S ld ) , (12) 

where V is the total system volume, S is the total interface area, and the geometric factor 7^ equals four times 
the ratio of the volume of a (d— l)-dimensional sphere to the surface area of a d-dimensional sphere of the same 
radius. On a d-dimensional hypercubic lattice of unit lattice constant, the number of bonds broken by the surface 
of a (i-dimensional sphere of radius R equals 2d times the corresponding discrete approximation to the volume of a 
(d — l)-dimensional sphere of the same radius. The interface area per unit volume, S/V, is therefore related to the 
total number of broken bonds, which is given by the sum in Eq. ( |10| ) and here called E, as S/V = 22/ (Nd'jd)- The 
relative error in this estimate comes from the discrete approximation to the (d— l)-dimensional volume and is ex 1/ R. 
The order-parameter dependent factor in Eq. (|l2j) is insignificantly different from unity for the times studied here, 
and it is therefore ignored. Inserting the expression for S/V in terms of E/7V in Eq. (113) and choosing d=3 yields the 
first term in Eq. (|l0|). The subtraction of unity is included to make R(t) vanish for the random initial condition. 
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FIG. 1. Average magnitude of the local order parameter at early times, ^{ip 2 (r,t)}, averaged over space and trials, shown 
vs time on a log-linear scale. The intial decrease at the earliest times is due to diffusional relaxation on short length scales. 
The near-exponential increase which follows is due to the local relaxation towards one or the other of the degenerate values, 
tp(r,t) ~ ±1. The solid curves are numerical solutions of the linear theory for the slightly anisotropic Laplacian used here, while 



the dashed curves correspond to the fully isotropic analytical result, Eq. 
with a very wide distribution of initial values, i[>o — 1. 
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FIG. 2. The characteristic length R from Eq. (10), shown vs time t, on a log-log scale. The dependence of R(t) on the 
system size indicates that finite-size effects influence the late-time growth exponent. A least-squares fit (the solid line) for the 
700 3 system, which should display minimal finite-size effects for the times shown here, yields an exponent estimate of n = 
0.511 ± 0.01 for t > 150. The dashed line is the predicted R(t) from Ref. ||, Eq. (0). 
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FIG. 3. Estimate of the instantaneous effective growth exponent vs the central time of the data used for the estimate. 
See text for details. Simple power-law growth with a constant exponent is seen only for t > 150 for 700 3 , the largest system 
considered here. Earlier times are clearly not associated with a constant exponent. 
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FIG. 4. A graphical test for t ' 2 growth in the late-time regime. The onset of finite-size effects occurs at progressively later 
times as the system size grows. For the 700 3 lattice the finite-size effects are unimportant for the times considered here. The 
straight line is a least-squares fit to the 700 3 data for t > 150, to the right of the large tick mark. 
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FIG. 5. The quantity A 1 (t) = 1 — -y/(^ 2 (r, t)) , which at late times is proportional to the characteristic length R(t). 
The solid line is a least-squares fit to the 700 3 data for t > 150, which yields an estimate for n of 0.512 ± 0.01. 



